Ladislas Nalborczyk & Thierry Phénix
Cours 1: Introduction à l'inférence bayésienne
Cours 2: Modèle beta-binomial
Cours 3 & 4: Modèle de régression linéaire
Cours 5: Markov Chain Monte Carlo
Cours 6: Modèle linéaire généralisé
Cours 7: Comparaison de modèles
Cours 8: Modèles multi-niveaux
Cours 9: Introduction à brms
Cours 10: Data Hackaton
On s'intéresse aux différences de taille entre hommes et femmes. On va mesurer 100 femmes et 100 hommes.
men <- rnorm(100, 175, 10) # 100 tailles d'hommes
women <- rnorm(100, 170, 10) # 100 tailles de femmes
t.test(men, women)
Welch Two Sample t-test
data: men and women
t = 2.1475, df = 197.71, p-value = 0.03297
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.2655634 6.2364837
sample estimates:
mean of x mean of y
173.0157 169.7646
On va simuler des t-valeurs issues de données générées selon l'hypothèse de non-différence entre hommes et femmes.
nSims <- 1e4 # number of simulations
t <- rep(NA, nSims) # initialising an empty vector
for (i in 1:nSims) {
men2 <- rnorm(100, 170, 10)
women2 <- rnorm(100, 170, 10)
t[i] <- t.test(men2, women2)$statistic
}
t <- replicate(nSims, t.test(rnorm(100, 170, 10), rnorm(100, 170, 10) )$statistic)
data.frame(t = t) %>%
ggplot(aes(x = t) ) +
geom_histogram() +
theme_bw(base_size = 20)
data.frame(x = c(-5, 5) ) %>%
ggplot(aes(x = x) ) +
stat_function(fun = dt, args = list(df = t.test(men, women)$parameter), size = 1.5) +
theme_bw(base_size = 20)
abs(qt(0.05 / 2, df = t.test(men, women)$parameter) ) # two-sided critical t-value
[1] 1.972035
alpha <- .05
abs(qt(alpha / 2, df = t.test(men, women)$parameter) ) # two-sided critical t-value
[1] 1.972035
tobs <- t.test(men, women)$statistic # observed t-value
tobs %>% as.numeric
[1] 2.147452
A p-value is simply a tail area (an integral), under the distribution of test statistics under the null hypothesis. It gives the probability of observing the data we observed (or more extreme data), given that the null hypothesis is true.
\[ p[\mathbf{t}(\mathbf{x}^{\text{rep}}|H_{0}) \geq t(x)] \]
t.test(men, women)$p.value
[1] 0.03297456
tvalue <- abs(t.test(men, women)$statistic)
df <- t.test(men, women)$parameter
2 * integrate(dt, tvalue, Inf, df = df)$value
[1] 0.03297456
2 * (1 - pt(abs(t.test(men, women)$statistic), t.test(men, women)$parameter) )
t
0.03297456
Comparaison de deux modèles:
\[ \underbrace{\dfrac{p(H_{0}|D)}{p(H_{1}|D)}}_{posterior\ odds} = \underbrace{\dfrac{p(D|H_{0})}{p(D|H_{1})}}_{Bayes\ factor} \times \underbrace{\dfrac{p(H_{0})}{p(H_{1})}}_{prior\ odds} \]
\[ \text{evidence}\ = p(D|H) = \int p(\theta|H) p(D|\theta,H) \text{d}\theta \]
L'évidence en faveur d'un modèle correspond à la probabilité marginale d'un modèle (le dénominateur du théorème de Bayes), c'est à dire à la likelihood moyennée sur le prior… Ce qui fait du Bayes Factor un ratio de likelihood pondéré (par le prior).
On lance une pièce 100 fois et on essaye d'estimer la probabilité \( \theta \) (le biais de la pièce) d'obtenir face. On compare deux modèles qui diffèrent par leur prior sur \( \theta \).
\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]
\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]
\[ \begin{align} \mathcal{M_{1}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(6, 10) \\ \end{align} \]
\[ \begin{align} \mathcal{M_{2}} : y_{i} &\sim \mathrm{Binomial}(n, \theta) \\ \theta &\sim \mathrm{Beta}(20, 12) \\ \end{align} \]
\[ \text{BF}_{12} = \dfrac{p(D|H_{1})}{p(D|H_{2})} = \dfrac{\int p(\theta|H_{1}) p(D|\theta,H_{1}) \text{d}\theta}{\int p(\theta|H_{2}) p(D|\theta,H_{2}) \text{d}\theta} = \dfrac{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(6, 10)\text{d}\theta}{\int \mathrm{Binomial}(n, \theta)\mathrm{Beta}(20, 12) \text{d}\theta} \]